Molecular docking, simulation and binding free energy analysis of small molecules as PfHT1 inhibitors

Antimalarial drug resistance has thrown a spanner in the works of malaria elimination. New drugs are required for ancillary support of existing malaria control efforts. Plasmodium falciparum requires host glucose for survival and proliferation. On this basis, P. falciparum hexose transporter 1 (PfHT1) protein involved in hexose permeation is considered a potential drug target. In this study, we tested the antimalarial activity of some compounds against PfHT1 using computational techniques. We performed high throughput virtual screening of 21,352 small-molecule compounds against PfHT1. The stability of the lead compound complexes was evaluated via molecular dynamics (MD) simulation for 100 nanoseconds. We also investigated the pharmacodynamic, pharmacokinetic and physiological characteristics of the compounds in accordance with Lipinksi rules for drug-likeness to bind and inhibit PfHT1. Molecular docking and free binding energy analyses were carried out using Molecular Mechanics with Generalized Born and Surface Area (MMGBSA) solvation to determine the selectivity of the hit compounds for PfHT1 over the human glucose transporter (hGLUT1) orthologue. Five important PfHT1 inhibitors were identified: Hyperoside (CID5281643); avicularin (CID5490064); sylibin (CID5213); harpagoside (CID5481542) and quercetagetin (CID5281680). The compounds formed intermolecular interaction with the binding pocket of the PfHT1 target via conserved amino acid residues (Val314, Gly183, Thr49, Asn52, Gly183, Ser315, Ser317, and Asn48). The MMGBSA analysis of the complexes yielded high free binding energies. Four (CID5281643, CID5490064, CID5213, and CID5481542) of the identified compounds were found to be stable within the PfHT1 binding pocket throughout the 100 nanoseconds simulation run time. The four compounds demonstrated higher affinity for PfHT1 than the human major glucose transporter (hGLUT1). This investigation demonstrates the inhibition potential of sylibin, hyperoside, harpagoside, and avicularin against PfHT1 receptor. Robust preclinical investigations are required to validate the chemotherapeutic properties of the identified compounds.


Introduction
Malaria is a major cause of morbidity and mortality despite concerted efforts to mitigate transmission [1]. While artemisinin-based treatment remains largely effective against malaria parasites, the potential emergence of multi-drug-resistant Plasmodium falciparum strains has necessitated the development of new therapeutic options [2][3][4]. To prevent malaria-associated public health crisis, there is an urgent requirement for novel antimalarial agents with high potency and favorable pharmacodynamic and pharmacokinetic profiles.
With advances in genomics and bioinformatics, it is reassuring that many new opportunities have emerged for the design and implementation of effective malaria mitigation strategies [5,6]. Computer-aided drug design (CADD) involves high-throughput screening of selective ligands to agonize or antagonize target structures [7,8]. CADD relies on the assumption that candidate compounds have affinity to protein targets with minimum side effects while having sufficient absorption, distribution, metabolism and excretion (ADME) properties [9,10]. Some promising malaria compounds that have recently progressed to clinical evaluation include KAE609 [11], M5717 [12], MMV390048 [13] and others [14].
P. falciparum-infected erythrocytes utilize up to 100 times more glucose than non-infected erythrocytes require because the parasite continuously metabolizes sugars from the host's erythrocytes to support its survival, growth, and replication [15]. The parasites have a significant survival advantage because of their transporter proteins which facilitate the capacity to convey a wide spectrum of sugar molecules successfully [16]. One of such transporter proteins is P falciparum hexose transporter (PfHT1) which is required for malaria parasite's survival and proliferation [16,17].
In this study, we tested the antimalarial activity of some compounds against PfHT1. We specifically performed structure-based high throughput screening of a library of 21,352 phytoligands following Lipinski's rules for potential small drug molecules [18]. We also performed molecular dynamics simulations on the lead compounds to elucidate protein motion by following their conformational changes through time. In all, we identified bioactive hit-to-lead compounds with varying potency and selectivity for PfHT1 over the orthologous human glucose transporter (hGLUT1).

Protein preparation
The 3-dimensional (3D) structure of Plasmodium falciparum Hexose Transporter 1 protein (PfHT1) was retrieved from Protein Data Bank (PDB) (https://www.rcsb.org). After considering residual factor (R)-value free, R-value work, R-value observed and overall resolution of the PfHT1 structures on PDB 6M20, 6ML2 and 6RW3, structure 6M20 was found to have the lowest values in all the parameters which indicated it could make a good target [19]. The PfHT1 protein with PDB ID 6M20 was refined by using Protein Preparation Wizard of Schrödinger-Maestro Release 2021-4 [19]. We assigned charges, bond orders and deleted water molecules to avoid inaccurately high binding scores [20]. Subsequently, hydrogens were added to the heavy atoms. The heavy atom root-mean-square deviation (RMSD) was fixed to 0.30Å using the optimized potentials for liquid simulations (OPLS) 2005 force field [21]. Lastly, we optimized amino acids using neutral pH. To determine the selectivity of our compounds as PfHT1 inhibitors versus, we downloaded the human GLUT1 (6THA) protein from the protein data bank https://www.rcsb.org, prepared the protein using the protein preparation wizard of Schrodinger suite. Subsequently, a receptor glide grid was generated and molecular docking of the receptor with our hit compounds was performed [22].

Ligand preparation
A total of 21,352 ligands from plants that have been documented to have antiplasmodial activities were used. The phytochemicals were downloaded in the structure-data file (SDF) format from the NCBI PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Ligand preparation was performed on the downloaded SDF files to assign proper bond orders and create a threedimensional geometry [23]. This was done in Maestro Schrodinger Suite 2017 using the Ligprep with OPLS 2005 force field [23]. In addition, ionization states were generated at pH 7.0 ± 2.0 with Epik 2.2 in Maestro Schrodinger Suite 2017 [24]. We generated 15 possible stereoisomers per ligand.

Receptor grid generation
The receptor grid was generated on the prepared protein. OPLS 2005 force field was used to generate the grids [25]. The van der Waal radii of the protein atoms were scaled by 1.0, the charge cutoff for polarity was 0.25. Furthermore, the receptor grid box was generated in each direction (x = 27Å, y = 27Å, and z = 27Å) and the box was set at the center of the cognate ligands with allowance for the binding pocket to accommodate any ligand [26]. The dock main after the grid generation was 27Å for each dimension (x, y and z).

Standard precision (SP) and extra precision (XP) ligand docking
Molecular docking experiment was performed by employing the glide executed in the Schrödinger suites [21]. The receptor was treated as a stiff structure while ligands were treated as flexible. The receptor grid was given a dimension suitable to accommodate ligand structures with a length � 14Å and a cubing docking grid was centered on Val318. The van der Waals scaling factor was set to 0.85 and 0.15 for non-polar atoms of the ligand and the partial charges limit value was set at -10.0 kcal/mole. High throughput virtual screening (HTVS) and standard precision (SP) scoring functions of glide were used and ligands were granted full flexibility. A post-docking minimization was carried out on output ligand-receptor complexes, reducing the initially collected 15 poses per ligand to five. The SP resultant compounds were further docked using extra precision (XP) mode with more accuracy and computational intensity. The configuration was without minimization, relaxation or flexibility. Based on the glide energy and XP glide rescoring, the procedure gave the lead ligand-receptor complexes. Subsequently, the glide module of the XP visualizer interface was used to examine specific interactions between ligands and proteins in addition to hydrogen bonds, hydrophobic interactions, internal energy, pi-pi (π-π) stacking interactions and RMSD.

Prime molecular mechanics with generalized born and surface area solvation (MMGBSA)
Binding free energy calculation was performed for the ligand-receptor complexes using the Schrodinger suite MMGBSA module integrated with Prime [21]. The binding free energy of XP Glide docked output complexes were evaluated using Prime MMGBSA. The evaluation of the complexes' relative energy was done with the OPLS3 force field and rotamer algorithm [27]. The free binding energy equation adopted was: ΔGbind = ΔGcomplex-(ΔGprotein + ΔGligand). A more negative score signifies a stronger binding energy [28].

Molecular dynamics simulation
To simulate the behavior of the biological environment, including water molecules and lipid membranes, we adopted molecular dynamics (MD) using Newton's equations to assess the motion of water, ions, tiny molecules, macromolecules or more complicated systems. To analyze the pattern of recognition of ligand-protein or protein-protein complexes, structural movements, such as those depending on temperature and solute/solvent, are crucial [29]. Molecular dynamics (MD) simulations were performed for 100 nanoseconds using Desmond Schrödinger [21,30]. Protein-ligand complexes used for molecular dynamics simulation were obtained from docking studies to provide a prediction of ligand binding status in static conditions.
Since docking is a static view of the binding pose of a molecule in the active site of the protein, MD simulation tends to compute the atom movements with time by integrating Newton's classical equation of motion [29]. The ligand binding status in the physiological environment was predicted using molecular dynamics simulations. Protein Preparation Wizard of Maestro Schrodinger Suite 2017 was used to preprocess the protein-ligand complex, which comprised complicated optimization and minimization [31,32]. All systems were prepared by the System Builder tool [33]. The OPLS_2005 force field was used in the simulation [34]. Solvent Model with an orthorhombic box was selected as transferable intermolecular interaction potential three points (TIP3P). We neutralized the models by adding counter ions where necessary. To mimic the natural physiological conditions, 0.15 M NaCl was added. Furthermore, the NpT ensemble with 300 K temperature and one atmospheric pressure was selected for complete simulation via Martyna-Tuckerman-Klein Barostat [35]. The models were relaxed before the simulation and the trajectories were saved after every 100 ns for analysis, after which the stability of simulations was evaluated by calculating the RMSD of the protein and ligand over time before analyzing the RMSF and protein-ligand contacts.

ADME-Tox properties
For the analysis of the pharmaceutical, physiological, biochemical, and molecular effects of the compounds, adsorption, distribution, metabolism, excretion, and toxicity (ADME-Tox) properties were calculated with the QikProp program of Maestro Schrodinger Suites [36]. The Qik-Prop predicted the physicochemical and pharmacokinetic properties of the compounds. It also assessed the tolerability of the analogues based on Lipinski's rule of five (that is, it does not violate more than one of the following criteria; no more than five hydrogen bond donors; no more than ten hydrogen bond acceptors; a molecular mass of less than 500 daltons; and a log P of less than five for octanol-water partition coefficient) [37].

ChEMBL validation of molecular docking
ChEMBL is an open bioactivity database containing binding, functional and ADMET information for a large number of drug-like compounds [38][39][40][41]. The bioactivity of PfHT1 was retrieved from CheMBL database (https://www.ebi.ac.uk/chembl/target_report_card/ CHEMBL4697/). Conical smiles containing 13,243 compounds were also downloaded from the database. The conical smiles data were viewed, cleaned and saved in Microsoft Excel 2016 as a comma-separated values (.csv) file. Using DataWarrior v.5.5.0 [42], the csv file was transformed to 2D (.sdf) format. Schrodinger 11.1 [31,36] was used to open and prepare the converted 2D (.sdf) file using ligprep (pH: 7, forcefield: OPLS3) [34]. The produced ligands were docked using the glide of target protein receptor with extra precision (XP) algorithm in Schrodinger 11.1. Subsequently, randomly chosen docking scores of 5000 compounds screened against PfHT1 in this study were plotted against corresponding inhibitory values obtained from the ChEMBL database after which the correlation coefficient was determined.

Data analysis
The raw trajectory files from the MD simulation run time and the Pearson correlation coefficient were generated in R (version 4.0.5) and visualized using "ggplot2" and "ggrepel" packages.

Molecular docking, MMGBSA/Prime binding energy of PfHT1-ligand complexes
We screened a library of 21,352 compounds against the target protein, PfHT1, and identified five hits from 437 compounds with excellent docking scores after thorough validation using the Lipinski rule of five (S1 File). To prove our hit compounds selected for PfHT1 over hGLUT1, we docked the five hit compounds into the binding pocket of hGLUT1 and found that hGLUT1 had very low affinity for the substrates. The molecular docking study of the five selected compounds and the target revealed that hyperoside had the highest glide docking score (-13.881 Å). The other four ligands; sylibin, avicularin, quercetagetin and harpagoside had -12.254 Kcal/mol, -11.952 Kcal/mol, -11.756 Kcal/mol and -11.258 Kcal/mol docking scores respectively (Fig 1). Furthermore, assessment of the ChEMBL-determined inhibition of PfHT1 and in-silico docking scores of PfHT1 indicated that the docking scores observed in this study were comparable with ChEMBL determined inhibitory values (Fig 2). For the residue interactions of a protein molecule with the ligand compounds we analyzed the protein-ligand complex structure and discovered that Val314 formed a hydrogen bond with four of the five complexes, while Ser317 and Gly183 appeared in three complexes with a hydrogen bond. Ser315, Asn48, Asn52, and Thr49 formed hydrogen bonds with two complexes each (Table 1). This suggests that the amino acid residues are essential for target binding.

Interaction of avicularin with PfHT1
Among the five hit compounds, avicularin showed a glide docking score that was closest to sylibin at -11.952 Kcal/mol. When the protein-ligand complex and the ligand atoms' contact with the target residues were observed, we found that the residue interactions of Ser315, Ser317, Asn48, and Val314 had double H-bonds (S5 Fig). The interaction involved back backbone and side-chain contacts, as well as hydrophobic contacts (Fig 4A).

Interaction of harpagoside with PfHT1
The PfHT1 target residues interacted with the atoms of the compound, the binding surface was controlled by a range of intermolecular interactions. The binding affinity depends on interactions at the bindings site and the non-specific forces outside the target binding region. The pattern of interaction between PfHT1 and harpagoside in the complex is shown in Fig 4B. The amino acid residue, Val314 formed triple pi-pi interaction with the ligand. We examined the interaction of harpagoside within the binding pocket of the target and discovered that the interaction was vigorous unlike hGLUT1 (Fig 5). This could be a result of the number of intermolecular interactions and the distance of the bonds. PfHT1 residues bound to the ligand through Asn48, Try49, Gly183, Ser317, and Asn52 made double H-bonds (S5 Fig). Asn52 formed H-bond back chain contacts with the harpagoside molecule ( Fig 4B).

Interaction of hyperoside with PfHT1
The interaction of hyperoside with PfHT1 is shown in Fig 1. The docking score of the PfHT1-hyperoside complex was -13.881 kcal/mol. The complex formed six conventional hydrogen bonds with Ser315, Ser317, Asn52, Gly183, and Val314 made a double hydrogen  bond with the ligand (Fig 4C; S5 Fig). The distance, 3 Å, encircled around the ligand was considered for screening. A pi-pi bond was also observed in Asn318 residue.

Interaction of quercetagetin with PfHT1
Quercetagetin had the fourth-highest glide docking score of -11.756 Kcal/mol and a good glide energy value. However, it is the only ligand with the least target residue contact. The targetligand complex interaction template showed it had double H-bond residue interactions with Val314 ( S5 Fig). We also observed pi-pi interaction via Asn318 and Asn48 amino acid residue. Meanwhile, the only pi-stack bond was on Ser317 (Fig 4D).

Interaction of sylibin with PfHT1
Compound Sylibin (CID5213) occupied the binding pocket of PfHT1 with the glide docking score of -12.254 Kcal/mol. Three hydrogen bond interactions were identified with the backbone amino acid residue Thr49, Gly183, and Val314 ( Fig 4E; S5 Fig). Gly183 and Val314 form H-bond contact with the backbone, while Thr49 forms an H-bond with the side chain. The interaction of the protein-ligand complex was robust as the ligand fits in perfectly into the binding pocket of the target.

Prime molecular mechanics with generalized born and surface area solvation (MMGBSA)
The prime MMGBSA integrated within the Prime Schrodinger suite was used to compute the binding free energy of the docked complexes. The relative free binding energies of sylibin, hyperoside, harpagoside, avicularin and quercetagetin were -75.43, -71.32, -63.62, -54.41 and -24.31, respectively as shown in Fig 1. The free binding energy further established the binding affinity of the selected ligands compared with the reference compound.

Molecular dynamics simulation
The selected ligands were evaluated for their conformational stability within the receptor's binding pocket. We examined the protein-ligand root mean square deviation (RMSD), protein root means square fluctuation (RMSF), ligand RMSF, protein secondary structure, proteinligand contacts and ligand torsion profile. The RMSD is the average deviation in the displacement of a group of atoms in relation to a reference frame for a given frame. Avicularin-receptor complex (lig-fit-prot) reached equilibrium after the first 30ns, the equilibrium was maintained till the end of the evolution with minimum and maximum values of 1.39 and 2.77Å, respectively. The Cα atoms reached a consistent fluctuation only after 0.5Å and was maintained all through the model (Fig 6). Likewise, the Cα atoms of sylibin, hyperoside, and harpagoside complex maintained an equilibrium state throughout the evolution time. There was a fair equilibrium in Cα atoms of the quercetagetin complex. Sylibin lig-fit-prot was observed to maintain a stable equilibrium for the period of 98 ns (Fig 6E; S6 Fig). While hyperoside and harposide lig-fit-prot reached equilibrium after 30 ns, the steady-state was maintained for 70ns. Quercetagetin complex, on other hand, showed a steady state from 15ns to 76ns, the oscillation was between 1.6 Å and 3.2 Å. We further observed a slight fluctuation between 76 and 90ns ( Fig 6D). Protein-ligand interactions were monitored throughout the simulation (Fig 7A-7E and S2  Fig). Hydrogen bonds play an essential role in protein-ligand binding. The conserved residues that formed hydrogen bond interactions were Asn48, Ser315, Lys51, Asn316, and Val 444. Interestingly, none of these residues interacted with quercetagetin, as shown in the RMSD and RMSF. Instead, quercetagetin showed remarkable RMSD and RSMF instability. Notable residues that form hydrophobic interaction with the ligands were Leu75, Leu81, Val443, and Val444, while more than 18 amino acid residues were conserved via water-bridge interactions with the ligands.

ADME-Tox evaluation
We evaluated the pharmacological and pharmacokinetic features of the hit compounds to predict their physiochemical characteristics. The properties represent the absorption, distribution, metabolism, and excretion (ADME) of the compounds. To estimate the ADME properties, we employed the Lipinski rule five (RoF: molecular weight (MW < 500); hydrogen bond acceptor (HBA < 10); hydrogen bond donor (HBD < 5); and predicted octanol/water partition coefficient (QPlogPo/w<5)) [37]. As shown in Table 2, sylibin, hyperoside and harpagoside did not violate any of the Lipinski rules. This makes the compound potentially druggable. About 75%, 62%, and 58% of hyperoside, sylibin, and harpagoside, respectively will be optimally absorbed into the system. The human oral absorption (HOA) of the three compounds showed cell permeability with considerable efficiency. Although avicularin had the highest HOA (79%), it violated one rule (hydrogen donor = 5). On the other hand, quercetagetin violated two rules (accptHB > 10 and donorHB > 5), and its HOA is very low compared with other compounds. This potentially suggests that the body can only absorb 10% of the quercetagetin.

Discussion
The obstruction of the glucose uptake pathway to starve out malaria parasites serves as a strategic way for new drug discovery. PfHT1 plays an essential role in the survival, proliferation and other metabolic activities of the parasite [15]. In this computational study, we screened a library of 21,352 compounds against PfHT1 and identified five phyto-ligand inhibitors of the protein. Interestingly, the five drug-like molecules identified have been previously reported to be active against some other diseases. For example, hyperoside has been reported to have neuroprotective [43], cardio-protective [44] and antioxidant activities [45,46]. Hyperoside (CID5281643), a quercetin3-O-D-galactoside (i.e., quercetin with a beta-D-galactosyl residue attached at position 3), is a flavonol glycoside present in a variety of vegetables and fruits [47]. It is predominant in Hypericum mysorense [48]. Avicularin (CID5490064) or quercetin glycoside, is a plant flavonoid with reported hepatoprotective property [49]. It has also been demonstrated to reduce C/EBP-activated GLUT4-mediated glucose uptake in adipocytes thereby inhibiting the formation of intracellular lipids [50]. It is isolated predominantly from Foeniculum vulgare and Juglans regia [49]. Harpagoside (CID5481542), on other hand, is a terpene glycoside chiefly from Harpagophytum procumbens (devil's claw). It has been reported to have anti-inflammatory properties against knee osteoarthritis [51]. Quercetagetin (CID5281680), a hexahydroxyflavone, is a plant metabolite that is predominant in marigold (Tagetes erecta) and Neurolaena lobata. It has reported antiviral [52]; antioxidant [53] and in vitro antilipemic potentials [53].
All five hits had similar interaction with the target compounds and fit into the target binding pocket with similar conformations, glide docking scores and very similar binding energies. The ligands had pi-pi interactions with Asn52, Asn318, Asn48, Ser317 and Phe53 residues and aromatic H-bond with Phe53, Asn48, and Val314 residues. Multiple H-bonds were present in all the docked complexes with residues Val314, Gly183, Thr49, Asn52, Gly183, Ser315, Ser317 and Asn48. Residues Val314, Ser317 and Gly183 formed an H-bond with 4 of the 5 ligandreceptor complexes. In addition, Asn48, Asn52, and Thr49 residues had H-bond intermolecular interaction with 3 of the 5 ligand-receptor complexes. These residues play important roles as they interact with ligands as hydrogen donor, hydrogen acceptor, and pi-pi interaction. Fonseca et al. [54] reported similar amino acids as essential residues in the binding pocket of PfHT1. The pi-pi interactions in most of the complexes were formed by Asn48, Asn52, Asn318, and Ser317. This observation is essential for drug development as most of the H-bond residues served as both donors and acceptors. The intermolecular features observed in this study could be explored to optimize ligand-receptor complexes. This can be utilized in the synthesis of entirely new molecules capable of interacting and inhibiting the target, with biological activity in vitro and in vivo [55].
MMGBSA analyzed the ligand-receptor intermolecular interactions by determining the ligand-receptor energy values and intermolecular interactions. With a high binding score and similar binding energy, the five compounds can bind PfHT1 receptors [56]. Our results have demonstrated a statistical correlation to experimental binding affinity when compared with extra precision glide docking score [57]. Taking together, the MMGBSA binding affinity, intermolecular pi-pi, and H-bond interactions with the conserved amino acid residues of PfHT1, sylibin, hyperoside, harpagoside, and avicularin are potential inhibitors.
MD simulations are crucial because they consider molecular structural motions which aid in the identification of hot spots, the interpretation of structural details at reported protein sites and the elimination of structural artifacts resulting from MD structural characterization conditions [29]. As a result of the robustness of MD simulation approaches, improved free energy estimates for protein ligand recognition can also be acquired and confirmed under experimental procedures [29]. To evaluate the stability of protein-ligand complexes for the hit compounds, 100 ns simulations were performed for each compound. Avicularin-receptor complex (lig-fit-prot) reached equilibrium after the first 30ns, the equilibrium was maintained till the end of the evolution with minimum and maximum values of 1.39 and 2.77Å, respectively. The Cα atoms reached a consistent fluctuation only after 0.5Å and was maintained all through the model. There was no drastic increase in RMSD value, possibly indicating stability of these two systems [58]. Quercetagetin complex, on the other hand, showed a steady state from 15ns to 78ns and the oscillation was between 1.6 Å and 3.2 Å. The slight fluctuation of the quercetagetin complex might be due to the number of rotatable bonds of the functional group which protrude outward of the target binding pocket [59]. Interestingly, there was no substantial fluctuation in the loop regions while comparing within models. The amino acid residues in avicularin, hyperoside, sylibin, and harpagoside model oscillated with little fluctuations; 0.5-1.5 Å, 0.5-1.6 Å, 0.5-1.0 Å, and 0.6-1.2 Å, respectively. We observed fluctuation in the quercetagetin model with 5.4 Å highest loop. The massive instability of the loops recorded was due to their inherent flexible nature of the ligand, which might be associated with the ligand interactions [60]. The instability of the quercetagetin model amino acid residues corroborates its lig-fit-prot RMSD, with large fluctuation between 76 and 90ns simulation run time. In essence, we found avicularin, hyperoside, sylibin, and harpagoside to be exceptionally stable in the binding active site of PfHT1 with negligible structural orientation and minimum conformational instabilities, rendering them likely inhibitors. Besides, avicularin, sylibin, hyperoside, and harpagoside showed excellent polypharmacological possibilities demonstrated by their docking scores, binding interactions, ADMET characteristics and interactions with receptor site residues. When we considered the protein-ligand interactions categorized into hydrogen-bonds, hydrophobic, ionic, and water bridges, we observed that all the selected complexes shared some amino acid residues that were conserved throughout the simulation run time. Moreover, the hit compounds were selective for PfHT1 over human orthologues (hGLUT1). This finding agrees with earlier reports of Joet et al. [61] which revealed that O-3 hexose derivatives inhibited uptake of glucose and fructose by PfHT1 when expressed in Xenopus oocytes. Selectivity of these derivatives for PfHT1 was confirmed by lack of inhibition of hexose transport by the major mammalian glucose and fructose transporters 1 and 5 [61]. The inhibition potential of compound 3361, an O-3 hexose derivative, has been effectively demonstrated in vitro, with high selectivity for Plasmodium spp. [62]. This potentially validates the compounds as drug candidates of interest.

Conclusions
Through molecular docking and MD simulation analyses, we have demonstrated that Asn48, Ser315, Ser317, and Val314 are likely essential amino acid residues for PfHT1 inhibition. We also showed that sylibin, hyperoside, harpagoside and avicularin are stable in the binding site and may efficiently inhibit PfHT1. Furthermore, the ADME properties of the four compounds make them potential druggable molecules which selectively inhibit the PfHT1 receptor over hGLUT1. These findings open a new line of investigation for in vivo modelling and preclinical assessment of the chemotherapeutic potentials of the identified compounds.